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Many modern production and measurement facilities incorporate multiphase systems at low pres¬ 
sures. In this region of flows at small, non-zero Knudsen- and low Mach numbers the classical 
mesoscopic Monte Carlo methods become increasingly numerically costly. To increase the numer¬ 
ical efficiency of simulations hybrid models are promising. In this contribution, we propose a novel 
efficient simulation approach for the simulation of two phase flows with a large concentration im¬ 
balance in a low pressure environment in the low intermediate Knudsen regime. Our hybrid model 
comprises a lattice-Boltzmann method corrected for the lower intermediate Kn regime proposed by 
Zhang et al. for the simulation of an ambient flow field. A coupled event-driven Monte-Carlo-style 
Boltzmann solver is employed to describe particles of a second species of low concentration. In order 
to evaluate the model, standard diffusivity and diffusion advection systems are considered. 


I. INTRODUCTION 

With technological advances in the 20th century, such 
as miniaturisation of micro electrical and bio-medical 
devices 0, 0, vacuum technology 0, space explora¬ 
tion 0 and production techniques requiring low pres¬ 
sure environments ii, flow systems have become of 
interest in which the mean free path A of individual fluid 
particles is not negligible in comparison to the geometry 
length-scale 1. Therefore, hydrodynamic models assum¬ 
ing a fluid continuum lose validity as non-equilibrium, 
kinetic effects gain importance. This transition is com¬ 
monly quantified by the Knudsen number Kn = X/l. In 
this publication, results of the development of a model 
for two phase flow systems with large concentration im¬ 
balances in complex geometries at low pressures are re¬ 
ported. These types of system are of importance e.g. 
for reduction of disturbances and quality improvement 
in the operation of advanced optical systems used in mi¬ 
crochip production 0 0 and chemical vapour deposition 
applications 0. 

While numerous numerical approaches capable of sim¬ 
ulating flows at finite Knudsen number, including mo¬ 
lecular dynamics ii , direct numerical simulation in¬ 
cluding slip corrections up to second order [10, discrete 
velocity methods m , lattice-Boltzmann methods [T^ . 
[T0 and Direct Simulation Monte Carlo (DSMC) 0 ex¬ 
ist, the latter is by far the most common and has been 
largely used to benchmark other approaches. Despite its 
success in modelling aerodynamics problems in the up¬ 
per atmosphere it has however become evident that for 
the case of flows in the lower transient and slip regime 
the DSMC approach still requires a very large number of 
particles to properly resolve a given flow problem [i0 . 

Being based on directly solving the Boltzmann form¬ 


alism as well, lattice-Boltzmann methods (LBM) too are 
in principle capable to describe gas flows in arbitrary 
flow regimes [ 10 . The commonly used D3Q19 single re¬ 
laxation time lattice-Boltzmann or lattice BGK model 
dealing with a linear approach to a discretised Maxwell- 
Boltzmann equilibrium distribution is however originally 
set up to recover the Navier-Stokes equations [ 10 . To 
extend its applicability to the simulation of low pressure 
gas flows in the slip flow regime of small non-negligible 
Knudsen numbers, boundary conditions can be extended 
to include slip at solid walls [I1,[IM0. For simulating 
flows up to moderate Knudsen numbers, phenomenolo¬ 
gical higher order lattice-Boltzmann models introducing 
wall function formalisms are used [l^, [2lj. Due to the 
existing high density gradients and dominating thermal 
velocity effects however, the LBM alone is not suited for 
the simulation of diffusive spreading of low concentrated 
components in rarefied gases. On the other hand the as¬ 
sumed low concentration of the second phase allows for 
a greatly simplified simulation approach, fully resolving 
particles where necessary, while maintaining computa¬ 
tional efficiency (detailed in Sec. IIIB|) . 

In order to increase the numerical efficiency of simu¬ 
lations at low Mach numbers/in the slip flow and trans¬ 
ition regime, several approaches involving the extension 
of existing models and (spatial) hybrid methods have 
been published. Fan et al. developed a modified DSMC 
method, called Information Preservation technique to re¬ 
duce statistical fluctuations associated with DSMC at low 
Mach numbers [ 10 . Chun and Koch reduced these fluc¬ 
tuations by solving the linearized Boltzmann equation 
together with a correction of the particle velocities after 
the collision step [10 . Hybrid approaches combining clas¬ 
sical Navier-Stokes solvers and DSMC are developed to 
reduce numerical costs [ 13 , [ 10 . In these approaches. 
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the computational domain is divided with respect to the 
local flow regimes and the DSMC method is only used 
in regions where non-equilibrium effects are dominant. 
Since the coupling between the grid-based Navier-Stokes 
solver and the particle based DSMC method is problem¬ 
atic, Burt et al. replaced the grid-based method with a 
simplified DSMC approach in the continuum region . 
Besides these hybrid approaches with increase in effi¬ 
ciency in mind other hybrid models have been developed 
to achieve a higher precision in approximations to solu¬ 
tions of the Boltzmann equation, such as e.g. MD-DSMC 
hybrid methods [13 • 

To extend the applicability of the hybrid approach to 
length and timescales immediately relevant to enginieer- 
ing problems regarding low pressure flows used to reduce 
contaminant concentrations, in the following a new hy¬ 
brid approach integrating LBM and DSMC simulation 
elements in a shared simulation domain is detailed. Here, 
the lattice-Boltzmann method as implemented in the ap¬ 
plication LB3D serves as an efficient (slip boundary aug¬ 
mented) Navier-Stokes solver for flows in complex geo¬ 
metries whereas a Monte Carlo particle model is used 
to study diffusion of the dilute contaminants throughout 
the system. The focus on contaminant transport in a low 
pressure environment allows to simplify the model signi¬ 
ficantly as here the global flow field can be considered 
without perturbation by the contaminants. Thereby the 
computational cost as compared to the hybrid models 
given above is further reduced- the contaminant mo¬ 
lecules are assumed to behave at equilibrium with the 
background-gas where the explicit simulation of particle 
representatives adds a thermodynamically consistent dif¬ 
fusion model to the athermal LB flow solver. 

After further elaborating on aspects of this newly 
developed approach, contaminant transport aspects in 
terms of thermal velocity, diffusivity and advection diffu¬ 
sion problems are evaluated. Results of simulations per¬ 
formed in a model of an experimental setup are reported 
before we close with concluding remarks. 


II. MODELLING ASPECTS 

The model aims to simulate dilute contaminants in a 
homogeneous atmosphere. This situation arises when low 
pressurised gas flows are employed to flush systems con¬ 
taining certain contaminants. The dilution is very im¬ 
portant for the simplification of the model to hold. In 
particular for such cases it is possible to neglect the action 
of the contaminant on the atmospheric gas. The LBM 
with intermediate Knudsen regime boundary corrections 
is employed to simulate the background gas. Single local 
contaminants are simulated as particles, where the colli¬ 
sion processes with the background gas are solved impli¬ 
citly by performing a Monte Carlo algorithm where the 
equilibrium velocity is Galileian shifted to account for 
advection. 

In effect, the simulation of the contaminant component 


is independent of the background gas model as long as 
the latter provides a (quasi-)static velocity field and ther¬ 
modynamic state as input. In the remainder we briefly 
revisit the modifications made to the LBM and give an 
overview of the Monte Carlo algorithm. Finally, the para- 
meterisation is discussed in the context of an example 
thermodynamic configuration. 


A. Lattice-Boltzmann methods for intermediate 
Knudsen flow 


The lattice-Boltzmann (LB) approach is based on the 
Boltzmann kinetic equation 




U • Vr 


/(r,u,t) = n, 


( 1 ) 


which describes the evolution of the single particle prob¬ 
ability density /(r, u, t), where r is the position, u the 
velocity, and t the time. The derivatives on the left hand 
side represent propagation of particles in phase space 
whereas the collision operator 17 takes into account mo¬ 
lecular collisions. 

In the LB method the time t, the position r, and the 
velocity u are discretised. In units of the lattice constant 
Ax and the time step At this leads to a discrete version 
ofEq. ID): 

fk{r + Ck,t+1) - fk{r,t) = flk, k = 0,l,...,B. (2) 


Our simulations are performed on a three dimensional 
lattice with R = 18 discrete velocities in direction of 
face and edge centres of a cube and a rest-vector of zero 
velocity per site (the so-called D3Q19 model). With a 
proper choice of the discretised collision operator 17 it 
can be shown that the flow behaviour follows the Navier- 
Stokes equation HI. We choose the Bhatnagar-Gross- 
Krook (BGK) form (H 

^k =-^(^fkir,t) - f^\u{r,t),p{r,t))^ , (3) 

which assumes a relaxation on a linear time-scale towards 
a discretised local Maxwell-Boltzmann distribution f^'^. 
Here, r is the mean collision time that determines the 
kinematic viscosity v = of the fluid. 

The equilibrium distribution function chosen 

within the scope of this work is element-wise given by 


/r = '^kP 


1 -b 4cfc ■ U-b i (Cfe • u)^ - 


1 

2 ^ 


Ul 


4 (Cfc • u) - i (Cfe • u) |u|‘ 


(4) 


herein is a weight-factor depending on the lattice geo¬ 
metry and Cfc are the lattice unit vectors and Cg is the 
lattice speed of sound. 

The particle velocity distribution / can be related to 
physical properties via its moments. Here follows e.g. 
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the macroscopic velocity vector u from the first order 
moment by 


u = - Vcfc/fc(x,t), (5) 

normalised by the local density 

B 

P = ( 6 ) 

k=0 


an accommodation parameter of a = 1 is chosen. This 
corresponds to diffuse deflection at a rough wall [l^ . 

As with higher Kn and/or increase in resolution the 
Knudsen layer occupies a range of several lattice sites 
depth, simple slip boundary conditions overestimate 
the velocity at the boundary. The reason for this is a 
significant reduction in the mean free path of particles 
in the vicinity of a surface, effectively lowering the 
local Knudsen number. Recently, several methods have 
been introduced to reflect this by the introduction of 
an effective mean free path Ag. Integrations with the 
lattice-Boltzmann method have been introduced by Tao 
et al., Hyodo et al. and Zhang et al. (2914^ . 


the zeroth order moment. 

As detailed in the introduction, for flows in systems 
either at very low pressure or of geometry on the nano¬ 
meter scale the characteristic dimensionless number, the 
Knudsen number Kn becomes none negligible. The mean 
free path for gases under ambient pressure (1.013 bar) is 
in the order of a few tens of nm, for classical technical 
flow problems in geometries much larger than this effects 
associated with higher Knudsen numbers can therefore 
be neglected. 

For Kn > 0.01 rarefication effects gain noticeable in¬ 
fluence on the flow characteristics and slip at the walls 
has to be taken into account [l^. This becomes relev¬ 
ant for macroscopically sized geometries in the regime of 
pressures of only a few Pascal considered here. Thus, to 
correct the lattice BGK (LBGK) Navier-Stokes solver for 
the discontinuous velocity field in the Knudsen layer, a 
slip boundary condition as proposed by Zhang et al. [l9| 
has been implemented. Other suitable slip boundary con¬ 
ditions were proposed by amongst others Tao et al. (^ . 
Ansumali and Karlin [l2l| and an extension of existing 
boundary conditions by Toschi and Succi introducing up¬ 
dated collision statistics by so-called virtual wall colli¬ 
sions [l^ . 

In the remainder of this work, a two-step extension of 
the boundary condition by a diffusive reflection regime as 
well as a wall function modifying the effective mean free 
path is used as proposed by Zhang et al.. Following their 
formulation, the Knudsen number in a LBGK model is 
given by 


Kn = 


[Yt-0.5 

V ^ I~ 


( 7 ) 


where L is the number of lattice sites used to resolve 
the characteristic length. From this it is immediately 
clear that the mean free path in this model is depend¬ 
ing on the relaxation time r only. It is thereby linked 
to the fluid viscosity. A first employed boundary condi¬ 
tion allows to control the reflective behaviour and thereby 
implicitly the slip by an accommodation parameter a 
defined between a = 0 relating to no-slip or bounce-back 
and 0 = 2 implementing full-slip or specular reflection. 
In the simulations described in the scope of this work. 


Here the correction proposed by Zhang et al. is ap¬ 
plied, formulating the effective mean free path 


l + 0.7e-‘^y/>'’ 

that contains a dependence on the distance from the 
nearest boundary node Ay. Using Eq. © and Eq. ([8]) 
the correction can be shown to correspond to a change in 
local kinematic viscosity, enterin g th e model via a (now 
local) relaxation time parameter [3ll j 


/Stt 

A 

V 8 

l-^0.7e-^y/^_ 


(9) 


B. A Monte Carlo Model for Contaminants 

To describe a second component present in very low 
concentrations only, a method derived from the DSMG 
approach is employed. As a solver to the Boltzmann 
equation, this model relies as well on the principal inde¬ 
pendence of the free movement and collision processes of 
particles. 

DSMG algorithms solve the Boltzmann equation in 
terms of free movement and collisions of particles rep¬ 
resentative of a thermodynamical state. Eree movement 
is calculated straightforward at a given (thermal) velo¬ 
city assigned to the particles. The distinction of the 
method lies in its approach to model particle collisions. 
After free flight for the duration of a given collision in¬ 
terval particles are binned and randomised collision part¬ 
ners in a bin volume are drawn. Each of these pairs 
is subsequently calculated to perform a momentum con¬ 
serving collision assigning a new particle velocity. This 
coarse grained collision model increases simulation effi¬ 
ciency dramatically as integration of the particle move¬ 
ment to predict explicit collisions is not necessary. At the 
same time the method has been shown to produce ther¬ 
modynamically accurate results, especially in the context 
of dilute systems with finite Knudsen numbers. 

In the approach presented here, the observation of col¬ 
lision time intervals is replaced by collision time calcula¬ 
tion from the mean free path travelled in the system. 







4 


Taking into account only atmosphere-contaminant in¬ 
teractions, atmospheric collision partners are not drawn 
from representative particles in a volume but rather cre¬ 
ated as pseudo-particles parameterised by properties of 
the LB held in the vicinity of a contaminant particle. 

An event driven algorithm focusing on collision events 
is implemented. Assuming local equilibrium, according 
to the mass of contaminant species C and temper¬ 
ature T the individual particle velocity is drawn from a 
Maxwell distribution. From the mean particle velocity 
the time span St to the next collision event is cal¬ 
culated as St = implicitly assuring the mean free 

path to be 




1 



( 10 ) 


In this derivation, properties of atmospheric species a 
were used. The particle density per unit volume Ua 
is calculated for an ideal gas. The symbols = 

{mam(^)/{nia + m(^) designate the reduced mass and ctqij 
the effective collision radius, calculated as the arithmetic 
average of atomic radii approximated by Lennard-Jones 
potential parameters given by Karniadakis et al. 0. 

At collision time the LBM lattice is used as binning 
grid. All contaminant particles present at a given lattice 
site undergo a collision with an ad hoc created pseudo 
particle^ generated as a representative of the local at¬ 
mospheric gas how. In particular, the atmosphere gas 
particles are assigned a velocity according to a Max¬ 
well distribution whose mean value u is shifted to rehect 
the velocity ulbm of the lattice-Boltzmann held. Per 
direction i this reads 




2TrkBT 


exp ■ 


—mg {uj — MLBM,i) 

2kBT 


( 11 ) 


Subsequent binary collisions between hard spheres are 
implemented classically, assuming conserved velocity of 
the centre of mass and thus global momentum conserva¬ 
tion, while updating the relative velocity of the collision 
partners . The centre of mass velocity and the rel¬ 
ative velocity of a contaminant particle C, and a pseudo 
particle a are given by 

u„ = -I-Uama)(m^-I-TO q)“\ (12) 


and 


- U„, (13) 

respectively. With u((, and u* denoting the post col¬ 
lision centre of mass and relative velocity respectively, 
kinetic energy and momentum is conserved in the col¬ 
lision by imposing = u™ and |u*| = |ur|. Since 
particles are assumed to be hard spheres all directions 
for u* are equally likely as scattering of hard spheres is 
isotropic. The post collision velocities of the particles can 
then be determined using a randomly chosen direction for 


the relative velocity considering |u*| = |ur| and the post 
collision equivalents of Eq. m and m- The updated 
contaminant velocity is employed for the calculation of 
the time interval St elapsed until the next collision event. 
The information of the pseudo particle is discarded. This 
collision process serves both to couple the contaminant 
particles to the atmospheric flow and as a thermostat (see 
section IIIIBI) . 

In the scope of this work the respective lattice- 
Boltzmann velocity field is pre-simulated until an equilib¬ 
rium state is reached. There exists however no principle 
limitation on the synchronisation of the relaxation pro¬ 
cesses. Rather the quality of the modelling of dynamics 
in the style of DSMC is improved by taking into account 
the local time-dependent equilibrium determined by the 
LBGK algorithm. Analytical approximations of trans¬ 
port coefficients in this regime suggest e.g. that the mu¬ 
tual diffusivity can under these assumptions be expressed 
as a Lorentz-approximation (mass-ratio) corrected self- 
diffusivity coefficient of the atmosphere gas (see section 

Imcl) . 

Rarefication effects are prominent in the Knudsen layer 
only. Thus special attention has to be paid to the bound¬ 
ary condition in the particle perspective as well. Since the 
main focus of the algorithm lies on true reproduction of 
the subscribed mean free path by timing collision events, 
the implementation of an accurate boundary condition is 
however straightforward. For the sake of simplicity again 
the lattice-Boltzmann site decomposition is used to define 
the boundary geometry, introducing discretisation errors 
in systems with curved boundaries. For these cases the 
method may be improved by another choice of surface 
dehnition, e.g. using interpolation techniques [12, 1^ . 
However, it has been shown that by resolving the Knud¬ 
sen layer by as little as four lattice sites the error ob¬ 
served in the flow field can be limited to the order of 5 
percent. For a resolution of 8 lattice sites the deviation 
reduces to some per cent. As discussed above for a work¬ 
ing model of flow in the intermediate Kn regime surface 
slip has to be taken into account as well as the reduc¬ 
tion in phase space volume in the vicinity of a bound¬ 
ary. Both requirements are met by a boundary condi¬ 
tion reducing a particle’s travelled distance if encounter¬ 
ing a wall. The algorithm evaluates, whether during free 
movement a particle is travelling into a different LB node 
and subsequently whether it is a boundary node. If the 
particle would enter a boundary node, its position is reset 
to the point of contact, effectively reducing its travelled 
path. The subsequent collision event is modelled as dif¬ 
fusive reflection inverting the velocity component normal 
to the boundary surface. The thermal velocity is drawn 
from the half set with an inverted velocity component 
normal to the wall. This treatment is implicitly redu¬ 
cing the mean free path proportional to the surface to 
volume ratio in the system as well as reducing the phase 
space volume. The boundary condition is evaluated for 
a simple advection diffusion problem for which an ana¬ 
lytical solution is available (see section UlIDI) . 
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C. Parameterisation 


The validity of the LBM is chiefly limited by the under¬ 
lying (thermo-) statistical principles and approximations 
in expansions. It is possible to define a LB specific Knud- 
sen limit based on the idea that the local equilibrium 
approximation on a lattice site does not hold anymore. 
Another way of putting this is given with the low Mach 
number limit which requires the transported momentum 
density on the lattice to be low. This clarifies that all 
modifications made to capture the intermediate Knudsen 
regime can only be phenomenological. As the evaluation 
of channel flow in the low Knudsen regime illustrates, 
modified continuum models can be suited to simulate 
flows in this regime. Nonetheless, when working with 
the modification of choice, leaving the core algorithm un¬ 
modified, the original restrictions have of course still to 
be observed. 


In order to be able to understand the units captured 
by the model, a conversion of the core units is instruct¬ 
ive. The scaling units of length Aa;, time At and mass 
Am are typically calculated by use of the constant speed 
of sound on the lattice as well as the kinematic shear vis¬ 
cosity imposed by the collision scheme. When aiming to 
model real systems, the usual choice of unit mass might 
be dismissed to calibrate the system in order to preserve 
numerical stability in cases were realistic pressures oth¬ 
erwise lead to very high densities or in cases where the 
relaxation rates are getting very low. 


In the cases employed here, starting out from length¬ 
scaling the lattice to match a physical system, length and 
time scales are fixed by comparing the speed of sound 
in the desired physical system with the one on the lat¬ 
tice. Keeping the lattice mass at unity, the overall mass 
scaling is then determined by comparison of the result¬ 
ing dynamic shear viscosity with the one of the desired 
physical system. 


A typical approach to the parameterisation of a sim¬ 
ulation system is made by deciding on the spatial resol¬ 
ution of a system, determining the lattice spacing Ax. 
Together with the lattice discretisation and the thereby 
determined speed of sound Cg this immediately fixes the 
time step At = Ax/cs and the conversion units for the 
kinematic viscosity For the case of a single fluid ad¬ 
hering to the ideal gas law p = pcg for the relation of 
pressure and fluid density, the consideration of the ther¬ 
modynamic state and modelled substance, i.e. temperat¬ 
ure, mass and pressure fixes a mass scale Am. Conversion 
of the dynamic viscosity rj then allows to determine the 
relaxation time or shear viscosity relaxation parameter 
T or Aj/ via = CgAt (-^ — i). Here, some tuning flex¬ 
ibility is given by the mass scaling, allowing to adjust 
the simulated mean fluid density against the relaxation 
time parameter. Table U gives some example numbers 
for arbitrary length scaling and Hydrogen at room tem¬ 
perature and 4 Pasc al p ressure obtained from the NIST 
chemistry webbook 341. 



Figure 1. Volumetric flow rate through a simple channel as a 
function of the Knudsen number, normalised by the flow rate 
measured at Kn=0.1. Reproduction of a validation published 
by Zhang et al. [3H| . The simulation results are obtained in a 
channel of 32 lattice units width. A finite slip boundary condi¬ 
tion was employed both alone and with a viscosity correction 
accountin g fo r a varying mean free path in the vicinity of a 
boundary [1^. HI. The reference values (symbols) are results 
of an exact solution to the linearized Boltzmann equation for 
hard spheres by Ohwada et al. [ssl. [3^. Using the combined 
boundary conditions good agreement can be obtained for the 
lower intermediate regime of Kn « 0.05..0.5. 


III. VALIDATION 


We start out with the re-evaluation of the implement¬ 
ation of the intermediate Kn corrected LBM boundary 
conditions as described by Zhang et al. [Ullli. Com¬ 
parison is made to results of the linearized Boltzmann 
equations in a narrow channel published by Ohwada et 
al. Non-dimensionalised values are compared. In 

a second part the thermal properties of the particle model 
as well as the coupling algorithm are tested. Through¬ 
out evaluations of the hybrid model SI units are used as 
they have been employed explicitly in the implementa¬ 
tion of the particle model. Here, we compare the simula¬ 
tion velocity distribution to the respective Maxwell dis¬ 
tribution both for the pseudo- and contaminant-particles 
over a range of temperatures. Third, the diffusive beha¬ 
viour of the model is checked. By focusing on a simple 
quasi-ID diffusion problem we ensure the validity of our 
choice of both mean free path A and diffusivity D. Fi¬ 
nally the lattice-Boltzmann velocity field, the coupling 
and the diffusivity model are integrated with our con¬ 
taminant boundary condition. Our results are compared 
to a solution of the transport equation in one dimension. 


A. Intermediate Knudsen numbers - The LBM 
implementation 

In order to assure the accuracy of the implementation 
of the boundary corrections by Zhang et al., flow at pre- 
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Property 

Formulation 

Physical example 

Length scale Ax 

Ax = x^'^^yx^^ 

1 m/1000 l.u. = 1.0 • 10"® m 

Time scale At 

At = c“/cf ■ Ax 

Ax/V3/1280S = 4.5105 ■ 10“'^ s 

Mass scale Am 

Am = {Ax^ ■ 

(Ax® ■ 2.445 • 10“® ^)/0.1 l.u. = 2.445 • lO”^"' kg 

Kinem. viscosity u 


3.581 = 1.6152 


Table I. Example numbers for arbitrary length scaling and Hydrogen at room temperature and 4 Pascal pressure obtained from 
the NIST chemistry webbook 


scribed Knudsen numbers in a simple channel is simu¬ 
lated. Evaluations are made both for slip boundary con¬ 
ditions only, as introduced in , and for combinations 
of a slip boundary condition and a viscosity correction 
to account for a locally varying Knudsen number at the 
boundary [MIlS- Simulations are executed in a pseudo- 
id channel of a size of 1x1x32 lattice units allowing for 
good numerical efficiency. The flow is driven by a body 
force ofF = l-10”^in lattice units. To assure to reach 
a steady state simulations are in all cases run for 100,000 
time steps. 

The flow profiles of the original publication are re¬ 
produced, reaching satisfactory agreement of simulated 
flow profiles and analytical solution to the linearized 
Boltzmann equation for hard spheres by Ohwada et 
al. [H, As depicted in figure |T] the combined bound¬ 
ary corrections allow to recover the theoretical flow rates 
within a few percent even up to low single digit Kn. 
The shown flow rate has been normalised by the re¬ 
spective flow rate assumed at Kn=0.1. Furthermore the 
Knudsen-Paradox is captured by the model. The min¬ 
imum flow rate measured in the simulations is reached 
around Knw 0.5 in good agreement with the results of 
the analytical calculations as well as other numerical res¬ 
ults reported by Toschi [I^ and Cercignani [s^. These 
measurements allow for confidence in the capability of 
the extended lattice-Boltzmann model to capture fluid 
behaviour in the intermediate regime. 



Figure 2. Probability density of contaminant particle velo¬ 
city for different system temperatures. Theoretical values are 
given by the Maxwell speed distribution for particles of a mass 
of 100 au at the respective temperature fEg. 1141) . We find ex¬ 
act agreement with the theory, verifying the correct operation 
of the coupling algorithm (Sec. Ill Bll . 


the collision algorithm, serving in addition to the coup¬ 
ling as a thermostat to the contaminants. 


C. The second law of Pick - Diffusion in a binary 
mixture with large density contrast 


B. The Maxwell speed distribution - The coupling 
algorithm 

To evaluate the coupling approach described in sec¬ 
tion IIIBI and ensure sound thermostatistic behaviour of 
the contaminant particles, the particle velocity statist¬ 
ics per collision event for a range of temperatures are 
measured. Parameters kept fixed are the ambient pres¬ 
sure p = 3 Pa, a contaminant mass of = 100 au and 
an atmosphere particle mass of ma = 2 au. Figure [2] 
documents the exact agreement between the simulation 
velocity field and the theoretical solution 

This result suggests correct thermalization of the pseudo¬ 
particle ’s velocity components as well as functionality of 


As detailed above, a main focus of application of the 
proposed model is the efficient simulation of diffusive 
transport in low pressure environments. The validity of 
the event driven algorithm, parameterised by the mean 
free path A is tested by comparison of the measured res¬ 
ulting diffusivity. The second law of Fick 

d 

—n^(a;) = Dacy'^riQ (15) 

defines the contaminant-atmosphere diffusivity in 
term of the second spatial and first time derivative of 
a concentration or density field. A solution in one di¬ 
mension gives the particle density n(^{x,t) at a locus x 
different from the initial position xq and time t in terms 
of the absolute particle number and Da<; as 


nc{x,t) 


Nc 


exp 


{x - I 


(16) 
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Figure 3. Particle density distribution at different times. All 
particles have been initialised in the point of origin. We find 
quantitative agreement with the theory (Pick’s second law, 
Eq. 1161) over a wide range of parameters. Here depicted is 
an example configuration; a system comprised of a Hydrogen 
atmosphere at p = 3 Pa and T = 295 K containing contam¬ 
inants of mass = 100 au. The diffusivity is calculated 
according to Eos. ll7lll8l to DqC ~ 1-31 • 10“^m^/s. 


FigureOillustrates the dynamics of the particle density 
field of a system comprised of a Hydrogen atmosphere at 
p — S Pa and T = 295 K containing contaminants of 
mass = 100 au. For the presented results, 1 million 
particles are initially placed at the origin. Using a mass 
corrected diffusivity 


Figure 4. Equilibrium particle density distribution in systems 
combining constant flow in positive x-direction and a diffus¬ 
ivity of DaC ~ 1.31 ■ Is. At xjxmax = 1 the system 

is delimited by a boundary acting on the contaminants (see 
section Hi Bll . The curves are given by the solution Eg. 1211 to 
the transport equation assuming an infinite particle reservoir 
at 1 > 1 and an open boundary at a: < 0. The simulation 
particle densities are re-normalised to the density measured 5 
lattice sites in front of the wall. Outside of the boundary layer 
we find the differential equation [19] excellently approximated. 
The fluctuating densities due to thermalised reflection at the 
wall do however not justify the infinite reservoir assumption 
of the theory. This results in strong variation of the obtained 
result. 

D. The transport equation - The particle boundary 
condition 
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3 / keT 
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with a mass correction factor 
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(18) 


excellent quantitative agreement of simulation and the¬ 
oretical predictions is found. The correction factor given 
by Eq. [TH| has allowed quantitative predictions for the 
diffusivity over a wide range of atmosphere/contaminant 
mass contrasts. It can be motivated by the non-trivial 
expectation values of momentum transfer at large mass 
contrasts. Even for simplified models of hard spheres, the 
asymmetry in momentum carried by species of different 
mass, travelling at different respective thermal speeds, 
requires higher order mass corrections. For a more in- 
depth discussion and exam ple calculations see the book 
of Chapman and Cowling [39|. The form of mass cor¬ 
rection reported here is again depending on the sim¬ 
plifying assumptions limiting interaction to atmosphere- 
contaminant collisions. For a mass contrast of one Cm 
reduces to unity. 


In this final evaluation of the model, the stationary 
state of a system described by the transport equation 

d 

—n(^{x) = —Dacy‘^nQ + Ua{x)VnQ=Q (19) 

is considered. Assuming the boundary conditions 

n^(0) = 0 ; nc_{x/xirLax = I) = 1, (20) 

the particle number density is normalised by a max¬ 
imum value assumed at a: = 1. This situation is de¬ 
scribing a system with a constant flow velocity in posit¬ 
ive x-direction and an infinite particle reservoir at x > I 
as well as an open boundary at x < 0. The infinite 
particle reservoir is here modelled by a boundary visible 
for the contaminant particles only. As before for the ex¬ 
ample system parameterisation a Hydrogen atmosphere 
at p = 3 Pa and T = 295 K containing contaminants 
of mass rrii^ = 100 au is selected. This corresponds to a 
diffusivity of DaC, ~ 1-31 • 10“^m^/s. Figure |4| depicts 
the normalised particle density over the normalised x- 
coordinate. Solid lines represent solutions to eg nation ITOl 
given by 



( 21 ) 














500 



Figure 5. Illustration of the flow path through the system. 
The velocities are not to scale. In the centre opening, the 
velocity is up to three orders of magnitude larger than in the 
remainder of the system. 


The simulation data is rescaled by the particle density 
5 lattice sites in front of the wall. This allows to com¬ 
pare the data to theory applicable to an infinite reser¬ 
voir. Measurements in closer vicinity to the wall suffer 
for the particle counts considered here from large dens¬ 
ity fluctuations introduced by the thermal wall boundary 
conditions. 


IV. SIMULATION OF CONTAMINANT 
SUPPRESSION BY BEZELS 

The model is applied to a setup where contaminant 
suppression by low speed flows through different openings 
at low pressure is measured. The simulations have been 
performed in a two-step process. First, a solution to 
the flow of the background gas was obtained by means 
of an adjusted LBM. Contaminants were subsequently 
simulated coupled to the quasi-static solution. In this 
part the LB parameterisation outlined in the preceding 
section HTCl is employed. 


A. System setup 

Figure [S] shows a ray-traced image of the simulation 
geometry used. The cylindrical shape of the chamber is 
reflecting an experimental setup. The system is resolved 
by 256 X 256 x 256 l.u.^, where in the centre-plate a round 
opening with a diameter of 14 mm is placed. 

The system geometry has been chosen to mimic a typ¬ 
ical experimental measurement setup. Here, openings are 
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Figure 6. Gas flow rate as a function of the pressure gradient. 
The LB simulation has been parameterised to the kinematic 
viscosity v = 3.581^ and speed of sound Ca = 1280 ^ of 
Hydrogen gas at a mean pressure of p = 4 Pa and temperature 
T = 295 K. A clear linear dependence of the flow rate Q on 
the pressure gradient is found even at very high flow speeds 
violating the low Mach number assumption (Ma « 0.27). 


introduced into the centre of a system with adjustable 
absolute pressure as well as relative pressure gradient. 
Flow rates are controlled by the pressure gradient and 
contaminants are injected on one side. Partial pressure 
is measured in both chambers. 

Single component fluid flow is parameterised by the 
choice of the mean lattice density p = 0.1 l.u. and a 
relaxation time of t = 5.3457 l.u. At the given discret¬ 
isation the viscosity wall-function (jH]) is found to have a 
non-zero value over 5 lattice sites. Full refractive bound¬ 
ary conditions are employed in the immediate vicinity 
of the bezel opening only. This has been established to 
enhance the stability of the LB in regimes of higher pres¬ 
sure gradients. This is justified as, while introducing an 
error in the exact form of the flow field throughout the 
system, flow rates in the central opening are not affected 
by this measure. 

The flow is driven establishing a flow rate by a Neu¬ 
mann condition [45| on the influx at a; = 0 against a 
fixed pressure Dirichlet condition [d^ at x = = 256. 

A snapshot illustration of the resulting flow field is depic¬ 
ted in figure [51 It is noteworthy that the flow rate in the 
region of the opening is up to three orders of magnitude 
higher than in the rest of the system, where vectors have 
been scaled equally to illustrate the flow path rather than 
the scales. 

In figure [5] atmosphere gas flow rates are plotted 
over the applied pressure gradients. The LB simula¬ 
tion has been parameterised to the kinematic viscosity 

= 3.581^ and speed of sound = 1280 f of Hydro¬ 
gen gas at a mean pressure of p = 4 Pa and temperature 
T = 295 K [ 3 ^ . The mean free path of hydrogen in this 
setup is calculated from tabulated data of the collisional 
cross section and mass to be Ac ~ 2.70 • 10 “^to, imply¬ 
ing a Knudsen number of Kn ~ 0.2 using the opening 
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diameter as limiting length scale. 

Linear dependence of the flow rate Q as measured in 
the volume of the opening on the pressure gradient is 
found. Even though in the simulations exceeding a pres¬ 
sure gradient of 1 Pa, the low Mach number limitation is 
clearly violated (at the highest value the Mach number in 
the peak flow is Ma 0.27), the functional dependence 
holds surprisingly well. This is due to the fact that this 
violation only occurs for a few lattice sites. The LB sim¬ 
ulations are run for 300,000 time steps after which the 
change in the field per time step is in the order of 10“^ 
and below. 



Distance [m] 


Figure 7. Snapshot of simulation data at the lowest and 
highest flow rates, respectively 100,000 collision events into 
the simulation where a total of 10,000,000 collision events 
was simulated. The data is a histogram of the amount 
of particles present in the respective lattice layer volume. 
The total number of particles was 100,000. The suppres¬ 
sion coefhcient is calculated from fitting a constant func¬ 
tion to these data in the ranges x £ {0.0050..0.0100} m and 
X £ {0.0175..0.0225} m. The step deviations in the central 
region reflect that the data is not normalised for the variation 
in local volume in the geometry. In particular, the opening 
in the region x £ {0.015..0.0152} m is clearly visible as min¬ 
imum. For the calculation of the suppression coefficient by 
Eq. (1231) particle numbers left and right of the opening are 
simply summed. 


B. Contaminant suppression 

The contaminant particles are initialised in the up¬ 
stream chamber in the range of a; = 0.0153 m to a; = 
0.0256 m. Each particle receives random spatial coordin¬ 
ates outside of the obstacle volume, as well as Maxwell- 
Boltzmann distributed random velocity components for 
a temperature of T = 295 K and particle mass of lOOart. 
The effective cross section is estimated from tabular val¬ 
ues for Hydrogen and Heptane (as an example of a heavy 


organic molecule of a weight of 100 a.u.) to be 0 

CTaC ~ I (2.91 • 10-1° -h 6.66 • 10-1°) m = 4.785-10-i°m, 

( 22 ) 

suggesting with equation OT a mean free path of the 
contaminant of approximately ~ 1.98 • 10~‘^Tn. Using 
the bezel opening diameter as typical scale, the obtained 
Knudsen number is Kn ^ 0.014. 

The contaminant dynamics are evaluated from concen¬ 
tration histograms of the system in the Cartesian co¬ 
ordinates. Eigure [7] illustrates these data in the main 
flow direction for the minimum and maximum flow rates 
employed {Q ~ 40 SCCM and Q w 500 SCCM, respect¬ 
ively) after 100,000 of 10,000,000 collision events. The 
difference in the system dynamics is already at this early 
stage distinctly visible. The concentrations and %- 
used to determine the suppression coefficient as 

5 = ^ (23) 

X 

are measured by fitting a constant in regions 
of undisturbed concentrations in the ranges x £ 
{0.0050..0.0100} m for the downstream and x £ 
{0.0175..0.0225} m for the upstream chamber. The devi¬ 
ations in the concentrations in the central region reflect 
that the data is not normalised for the variation in local 
volume in the geometry. In particular, the opening in 
the region x £ {0.0150..0.0152} m is clearly visible as 
minimum. 

The suppression coefficients obtained for the various 
flow rates are presented in figure [71 In agreement with 
findings involving an earlier model, an exponential de¬ 
pendence of the suppression on the flow rate is found [13 . 

This test case illustrates the applicability of our hy¬ 
brid model to real world systems with industrial relev¬ 
ance. Quantitative comparison with experimental data 
is a natural next step. 

V. CONCLUSION AND OUTLOOK 

Motivated by industrial applications requiring contam¬ 
inant suppression in low vacuum environments, a hy¬ 
brid simulation model has been refined and validated. A 
Monte Carlo algorithm to introduce a diffusivity model 
to a LBM setup has been developed and improved in 
thermostatistic context. 

The problem considered comprises an ideal gas which 
is flushing a system pressurised in the single digit Pas¬ 
cal range, where medium sized organic contaminants are 
present in the system. The sparsity of the contaminant 
molecules allows complete neglect of their impact on the 
gas flow. As a result it is necessary to couple the solution 
of the gas flow into the simulation of the contaminants 
only. 

This important simplification allows efficient calcula¬ 
tion of the diffusive transport (The contaminant simu¬ 
lations are executed as a single process and finish in a 
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Figure 8. Measured suppression coefficients over varying flow 
rate through the opening. The suppression behaviour of the 
system can be approximated by an exponential relation to the 
flow rate observed in the system. 


matter of hours). By only requiring a quasi static solu¬ 
tion to the flow field, the Monte Carlo portion of the 
model is furthermore independent of the choice of flow 
solver used to obtain said solutions. 

With the intermediate goal of simulating and predict¬ 
ing flow on scales relevant to engineering, the LBM has 
been selected due to its adaptability to complex geomet¬ 
ries. While the method is not readily applicable to in¬ 
termediate Knudsen flow regimes, phenomenological cor¬ 
rections can be introduced to recover the flow behaviour 
observed in rarefied systems with sufficient accuracy. 

The Monte Carlo model is distinguished by combina¬ 
tion of two central ideas. First, the coupling of a local 
advective velocity by shift of the expectation values of 
the Gaussian thermal velocity distributions constituting 
the local Maxwellian. Second, the parameterisation of a 
collision event-driven algorithm determining time scales 
by the mean free path and equilibrium velocities. 

The resulting algorithm has been validated to recover 
correct thermal velocity distributions as well as super¬ 
posed advective velocities. Since the coupling algorithm 
implicitly acts as a thermostat, the Maxwell-Boltzmann 
distribution of particle velocities is assured at all times. 

Using collision properties of hard spheres with an ad¬ 
ditional mass contrast correction, quantitative prediction 
of the diffusivity of the contaminants for a given thermo¬ 
dynamical state and particle mass has been successful. 
The mass correction term found here may be of prac¬ 


tical relevance to more than this application as it illus¬ 
trates and quantifies the relevance of mass contrasts in 
momentum transfer and resulting diffusivity. 

A more complex test case combines an effective 
advection-diffusion, or transport model with a test of the 
boundary conditions. It has to be refrained from devel¬ 
oping a quantitative solution to the particle distribution 
in the boundary layer. Quantitative agreement with the 
approximation of an infinite fixed concentration source 
could however be obtained when a density integrated over 
the boundary layer was used. 

In a final step the developed combined method has 
been applied to suppression of contaminant diffusive 
transport by gas flows through an opening in the low 
intermediate Knudsen regime. Qualitative evaluation of 
gas flow rate and resulting suppression has yielded prom¬ 
ising results. As expected there exists a well defined lin¬ 
ear dependence of the flow rate on the applied pressure 
gradient. Furthermore the resulting suppression is an 
exponential function of the flow rate. 

Further development and application of the model can 
be envisioned in both experimental as well as theoretical 
context. True to the original motivation of the work a 
natural next step is the comparison with experimental 
results. To this end, this work provides a parameterisa¬ 
tion framework both of the LBM and the Monte Carlo 
algorithm where input can be provided in SI units. While 
some limitations such as the low Mach number limit and 
the principal continuity assumption apply, the method is 
able to cover a wide range of technically relevant scales. 

Theoretical extensions of the model can be made by 
means of refining the algorithm to increase its overall 
accuracy. A very obvious flaw of the current model liv¬ 
ing on the LBM lattice is the very limited resolution of 
boundary shapes. Also the coupling to models capable 
of providing exact solutions to the Knudsen layer flow to 
some problems could be very interesting. Also extension 
of the collision model to include soft potentials, etc. is of 
interest. 
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